Stability and metastability of clusters in a reactive atmosphere: 
theoretical evidence for unexpected stoichiometrics of ISAgnjO^ 
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By applying a genetic algorithm in a cascade approach of increasing accuracy, we calculate the 
composition and structure of MgMOa; clusters at realistic temperatures and oxygen pressures. The 
stable and metastable systems are identified by ah initio atomistic thermodynamics. We find that 
small clusters (M < 5) are in thermodynamic equilibrium when x > M . The non-stoichiometric 
clusters exhibit peculiar magnetic behavior, suggesting the possibility of tuning magnetic proper- 
ties by changing environmental pressure and temperature conditions. Furthermore, we show that 
density-functional theory (DFT) with a hybrid exchange-correlation (xc) functional is needed for 
predicting accurate phase diagrams of metal-oxide clusters. Neither a (sophisticated) force field nor 
DFT with (semi)local xc functionals are sufficient for even a qualitative prediction. 

Keywords: Clusters, Ab initio Atomistic Thermodynamics, Genetic Algorithm, DFT, Magnesium Oxide, 
Replica-exchange Molecular Dynamics 



In the search for novel functional materials, atomic 
(sub)nanometer clusters are widely studied as model sys- 
tems exhibiting unique size-dependent properties often 
even qualitatively different from bulk materials. For ex- 
ample, small clusters may exhibit completely new lo- 
cal structures, stoichiometrics, electronic and magnetic 
properties unknown in the bulk materials [T]. Heteroge- 
neous catalysis is just one important example where all 
mentioned issues are of fundamental importance [JHH] . 

The composition and structure of clusters are deter- 
mined by thermodynamics and kinetics at the relevant 
temperature (T) and the nature of the environment. In 
thermodynamic equilibrium, only structures and compo- 
sitions that minimize the free energy of the combined 
gas-l-cluster system will be stable. Although a system 
is often not in thermodynamic equilibrium, thermody- 
namic phase diagrams serve as guidelines and important 
limits for predicting properties and functions of real ma- 
terials. In this Letter, we address the issue of stability 
and metastability using a model system that is relevant 
for many practical applications: free metal (Mg) clusters 
in an oxygen atmosphere. 

Most of the previous research on clusters focused on 
properties of stoichiometric (MgO)j\/ clusters [3 [5HT5]. 
and only few attempts have been made to study prop- 
erties of the non-stoichiometric [Mg^Oa,] clusters. [BHTU] 
However, the decisive issue of stability and metastability 
of clusters with different compositions at realistic condi- 
tions (exchange of atoms with an environment) has not 
been addressed so far. 

We consider a wide range of MgA/Oj; cluster sizes: 
1 < M < 10 and x determined by thermal equilibrium 
with the environment at given temperature T and par- 
tial oxygen pressure P02 • For each stoichiometry, the 
energy is minimized with respect to both geometry and 
spin state. Very unexpectedly, our results reveal (see Fig. 
[1]) that non-stoichiometric clusters with x > M are more 
stable at realistic (T, P02) when M < 5. For bigger 



clusters the stoichiometric composition, (MgO)A/, is pre- 
ferred, although with geometries very different from the 
bulk rock-salt structure. 

For reliable determination of low-energy structures, 
we have developed a massively-parallel genetic algorithm 
(GA). GA mimics the process of "natural" selection to 
evolve a set of atomic structures until the structures that 
fit best chosen selection criteria are found [IBl - fIS] . In 
this work, we search for structures that minimize the 
DFT total energy within each stoichiometry. Details of 
the implementation and validation of our GA scheme are 
given in Ref. |19| . It is based on an initial pre-screening 
of structures obtained by a GA performed with a (re- 
active) force field (FF-GA). The approach proceeds in 
terms of a cascade in which successive GA runs employ 
higher levels of theory, with each next level using infor- 
mation obtained at the lower level. We also note that the 
method is massively-parallel: new structures are gener- 
ated by combining two structures randomly selected from 
the running pool, via a set of modified cut-and-splice 
operators [121 HO] by many parallel processes (replicas) , 
each replica optimizes its newly generated structure, and 
a pool of structures shared among all the replicas is 
updated each time a replica finds a new optimized struc- 
ture. With the communication limited only to writing 
to / reading from the shared pool of structures, one can 
then run as many parallel replicas as CPU availability 
permits. 

The free energy as a function of T and P02 is calculated 
for the minimum of the potential-energy surface (global 
minimum, GM) and (energetically) adjacent local energy 
minima for each stoichiometry using the ab initio atom- 
istic thermodynamics (aJAT) approach pil - Bl] . recently 
extended to cluster systems j25j. The thermodynamic 
phase diagram is then constructed by selecting cluster 
compositions and structures with the lowest free energy 
as a function of (T, pQ^ ) . 

The set of highest binding-energy structures (within 
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FIG. 1: Free energy of formation of thermodynamically most stable non-stoichiometric (MgA/O^ with M ^ x) relative to 
stoichiometric (M — x, denoted as (MgO)A/) clusters at po^ — 1 atm and various temperatures. The geometries were 
optimized with PBE+vdW, and the electronic energy was calculated using PBEO+vdW. The global-minimum structures 
and spin multiplicity A4 (indicated as exponents to the chemical formula of the clusters, if different spin states coexist) at 
T = 300 K are also shown. Structures with chemical formulas in green are the most stable at that size M at T = 300 K, but 
other stoichiometrics coexist within 2k^T (i.e., the minimum concentration of clusters of a given competing type is at least 
10% of the total mixture). At other temperatures, the most stable stoichiometry and/or isomer may differ from the one at 
T = 300 K. The structure labeled in magenta, (MgO)io, is present as three coexisting isomers (the one shown is the most 
stable, more details can be found in Ref. I19p . Coordinates for all structures are available in the Suppl. Material. 



4 eV from the GM) obtained from FF-GA is used as 
initial pool for a DFT-based GA structure search. In 
this way, the initial pool contains ^ 30 structures for a 
small system like Mg202 and ^ 300 for the largest sys- 
tem we report here, MgioOi5. We label this linked ap- 
proach as FF/DFT-GA. Starting from FF-relaxed struc- 
tures, rather than from initial random structures, greatly 
reduces the overall computational time of the procedure 
(up to a factor of ten for our systems). This is because 
the geometry optimization of structures at the DFT level 
is computationally the most expensive part of an ah initio 
GA. A FF-GA search risks however to introduce a bias 
that would affect the results obtained via FF/DFT-GA. 
We checked T^ that this is not the case for the chosen 
force field (reaxFF [5^1 [57]), by verifying that (expen- 
sive) DFT-GA runs, i.e., started from randomly gener- 
ated pools, yielded exactly the same final list of isomers 
obtained through a (much less expensive) FF/DFT-GA. 
It must be stressed that just locally relaxing via DFT the 
structures found by the FF-GA is insufficient for an ac- 
curate sampling of the DFT PES, especially at large O2 
coverages. In fact, we observed that by adopting this lat- 
ter procedure one can miss some low-energy isomers |19j . 

A challenging problem of any GA is to guarantee that 
the lowest-energy structure found by the algorithm is 
indeed the GM. We address this problem by applying 
replica-exchange molecular dynamics (REMD) [25l f28l 
129] , a (computationally very expensive) reference method 
that explores the PES using the dynamics of the system 
at different temperatures simultaneously, and is exhaus- 
tive if the system is ergodic. Our implementation of GA 
gives indeed the same GM as REMD [19] . 

Within DFT-GA, the actual "cascade" unfolds at the 
local optimization step, which is initially performed with 
lower-level (computationally relatively cheap) numerical 
settings. Next, structures that are not previously found 



(details on the comparison of structures are given in 
Ref- IE]) and have energies within 2.5 eV from the 
current GM candidate, are further relaxed using higher- 
level settings [3D|. All DFT calculations have been per- 
formed with the FHI-aims package [31]. For the lower- 
level energy minimization, the van-der-Waals (vdW)- 
corrected [31j PBE [33] xc functional (PBE+vdW) and 
"light" numerical settings with (numeric atom-centered) 
basis-functions set "tier 1" were used [31 , and forces were 
converged to better than 10~'^ eV/A. At the higher level, 
energy minimization is performed with PBE-f vdW func- 
tional, "tight - tier 2" settings, and forces were converged 
to better than 10^^ eV/A. The energy of these further op- 
timized structures are evaluated via vdW-corrected [35] 
PBEG ^ hybrid xc functional (PBEO+vdW), with 
"tight - tier 2" settings [35]. Within DFT-GA, the latter 
energy is used to evaluate the fitness function, i.e., [19] 
a mapping of the energy interval between highest- and 
lowest-energy cluster in the running pool into the inter- 
val [0, 1]. The higher the value of the fitness function for 
a cluster, the higher is the probability of selecting it for 
generating a new structure. 

It has been recently shown "35" that HSE06 func- 
tional [32; yields a good description of the level align- 
ment and electron transfer in MgO. Based on compari- 
son for selected clusters, we find that PBEG, which be- 
longs to the same family of functionals, yields forma- 
tion energies essentially identical to IISE06 (within 0.05 
eV for MgA/O-c neutral clusters, see also jTO]) The ac- 
curacy of PBEO+vdW energies for MgMO^; clusters was 
also validated against the highest level currently achiev- 
able within the DFT framework, i.e., the renormalized 
second-order perturbation theory (rPT2) [38j . here ap- 
plied on PBE orbitals (rPT2@PBE) [35]. 

We find that PBE+vdW strongly overestimates sta- 
bility of clusters with larger x, resulting in a qualita- 
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FIG. 2: Energy of 02-adsorption on MgOo; clusters (energy 
of the reaction MgOi, + O2 — >■ MgOa;+2 calculated at the 
PBE+vdW GM geometry) calculated with different xc func- 
tionals and reaxFF 26 . For the cases where singlet (S) and 
triplet (T) spin states are almost degenerate, the two sets of 
energies are shown (for DFT only, since ReaxFF does not 
describe spin). In these cases, the spin state of MgO^: is cho- 
sen to be singlet for triplet Mg02,+2, and triplet for singlet 
MgOi:+2, to preserve the spin-conservation rule (O2 is always 
kept in its ground triplet electronic state). The geometry dif- 
ferences of clusters with different spin states are invisible at 
this scale. 



tively incorrect prediction that O2 adsorption would be 
favored over desorption up to a large excess of oxy- 
gen. Such behavior is not confirmed by hybrid func- 
tional or rPT2@PBE. To demonstrate this, the DFT 
energy of the reaction MgO^; -I- O2 — J' MgOj:+2 calcu- 
lated with different functionals is shown in Fig. [2J The 
PBE+vdW overestimation of energies at large coverages 
implies that, during the GA scanning, a fitness function 
based on PBE+vdW energy would lead to a biased sam- 
pling, because low-energy structures are different at the 
converged PBEO+vdW level. As can be seen in Fig. [2] for 
lower 02-coverage the difference between PBE+vdW and 
PBE0+vdW/rPT2@PBE energies of O2 adsorption on 
Mg and Mg02 is small despite the error in the O2 binding 
energy (the calculated O2 binding energy is 6.23 eV for 
PBE+vdW, 5.36 for PBEO+vdW, 4.59 for rPT2@PBE, 
and the experimental value is 5.21 [4^). This can be ex- 
plained by error cancellation for the clusters: If an O2 
molecule adsorbs non-dissociatively on MgMOa;, the er- 
ror in the description of MgM0a;+2 will cancel with the 
error in the description of O2 when calculating the ad- 
sorption energy. Indeed, we find that adsorption of O2 on 
Mg and Mg02 does not lead to breaking of 0-0 bonds. 
The difference between PBE+vdW and PBEO+vdW en- 
ergies of O2 adsorption on MgO for 7M = 3 is due to 
the difference in the description of the singlet state of 
MgO itself. For clusters with a; > 5, correction of the 
O2 binding energy error increases the difference between 
PBE+vdW and PBE0+vdW/rPT2@PBE adsorption en- 
ergies (see the Suppl. Material) . The tendency of PBE 
(and, as a consequence, of PBE+vdW) to overbind O2 
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FIG. 3: Phase diagram for Mg40a; clusters in an oxygen at- 
mosphere. The geometries are optimized with PBE+vdW 
and the total energies are calculated with PBEO+vdW. The 
sand-colored unlabeled areas are regions where different com- 
positions (at least the adjacent ones) coexist (see text). The 
thick black line is the 0-rich limit HH EH [21]: Above this 
line, O2 droplets condensate on the clusters. The rectangle 
encompasses a region accessible to experiments. 



molecules at high coverage holds also at larger M (see 
the Suppl. Material). 

At given T, P02: ^-nd M, the stable stoichiometry 
of a MgMO:r cluster is determined via azAT, i.e., by 
minimizing the Gibbs free energy of formation [191 I41j : 
AGfiT,poJ = Fms^oAT) - Fus^iT) - x,jLo{T,po,). 
Here, -FMgj,jO^(r) and Fyig,^iT) are the Helmholtz free 
energies of the Mga^Oa; and the pristine MgM cluster 
(at their ground state with respect to geometry and 
spin), respectively, and pi.oiT,po^) is the chemical po- 
tential of oxygen. As explained in [HI [41], Fyig^jO^iT) 
and -pMgjj- (T) are calculated using DFT information and 
are expressed as the sum of DFT total energy, DFT vi- 
brational free energy in the harmonic approximation, as 
well as translational, rotational, symmetry- and spin- 
degeneracy free-energy contributions. The dependence 
of /io on T and po^ is calculated using the ideal (di- 
atomic) gas approximation with the same DFT func- 
tional as for the clusters [TH]. The phase diagram for 
a particular M is constructed by identifying the lowest 
free-energy structures at each (T, ^02)- As a represen- 
tative example, we show in Fig. [3] the phase diagram for 
M = 4. Phase diagrams based on reaxFF, PBE+vdW, 
PBEO (without vdW), and rPT2@PBE can be found 
in the Suppl. Material. We find that at all DFT lev- 
els, the phase diagrams are qualitatively and quantita- 
tively very similar for T > 200 K and po^ < 10~^ atm. 
For higher pressures and/or lower temperatures, how- 
ever, PBE+vdW predicts larger x in MgjyiOx as thermo- 
dynamically more stable, compared to PBEO+vdW and 
rPT2(8)PBE. This is consistent with the results shown in 
Fig. [2] i.e., PBE tends to favor adsorption of a larger 
number of O2 molecules. PBEO and PBEO+vdW yield 
almost identical phase diagrams. Thus, at the consid- 
ered cluster sizes the vdW interactions within the scheme 



of Ref. [52] do not play a crucial role. ReaxFF-based 
diagrams are very different from DFT-based ones and 
no conclusions, even qualitative, can be drawn from the 
analysis of the reactive force field results. Furthermore, 
we find [lOj that neglecting translational, rotational, and 
vibrational, symmetry- and spin-degeneracy-related free- 
energy contributions results only in slight changes in the 
phase diagrams, comparable to the differences between 
PBEO-fvdW and rPT2@PBE diagrams, for aU consid- 
ered M. 

We have constructed phase diagrams for 1 < M < 10. 
At thermodynamic equilibrium in a wide range around 
normal conditions, small Mgjv/O^: clusters (M < 5) are 
found to be preferentially non-stoichiometric [x > M). 
For sizes Af = 6, 7 we observe a competition between 
stoichiometric and non-stoichiometric stable structures, 
while for sizes M > 8 the stoichiometric composition is 
the most stable, with the energy gap between the stoi- 
chiometric and non-stoichiometric structures increasing 
with size M . In Fig. [T} we also report the relative sta- 
bility of stoichiometric and non-stoichiometric structures 
at low (T = 100 K) and high (T = 1 000 K) tempera- 
tures. The relative stability of stoichiometric structures 
increases with temperature. We also find (see Suppl. Ma- 
terial) that lowering po^ shifts the equilibrium towards 
stoichiometric structures. 

In Fig. [1} we report also the spin multiplicity M of 
the low- free-energy structures. We find that in general 
MgMOa; non-stoichiometric clusters have several near- 
degenerate (within ^0.05 eV) electronic states with dif- 
ferent spin. The highest multiplicity we find is A^ = 7 for 
the case Mg40i2, and values of A^ = 3 and 5 are largely 
represented. The analysis of the spin-density shows that 
the unpaired electron density is localized on O2 and O3 
moieties (see examples in the Suppl. Material). The or- 
bitals occupied by the different unpaired electrons in the 
same cluster have vanishing overlaps due to large sepa- 
ration. The different spin-states for each isomer can be 
formed and coexist in oxygen atmosphere without vio- 
lating spin conservation rules, since successive adsorp- 
tion and desorption of the (triplet) O2 molecules in the 
gas phase allows the clusters to reach all the stable spin- 
multiplicities observed for each cluster. The thermody- 
namically favored access of oxygen in small clusters is in 
sharp contrast to bulk pristine MgO, where stoichiomet- 
ric composition is strongly favored. In all these clusters, 
the coordination of O atoms to Mg is maximum two. 
When the number of Mg is increased such that every 
O atom can be coordinated to three or more Mg, and 
at the same time the clusters have large HOMO-LUMO 
gap (i.e., all valence Mg electron are transferred to O^;), 
then stoichiometric clusters become thermodynamically 
stabilized. 

All the quasi-degenerate spin states are populated at fi- 
nite temperature, thus all the non-stoichiometric clusters 
with energetically quasi-degenerate states are paramag- 



netic. In contrast to non-stoichiometric clusters, we find 
that stoichiometric structures are always singlet, sepa- 
rated from the higher-multiplicity states by at least 1 eV. 
Thus, the stoichiometric clusters are diamagnetic. By 
looking at Fig. |3j we note that in the range of real- 
istically achievable pressures encompassed by the rect- 
angle, non-stoichiometric (paramagnetic) structures are 
thermodynamically more stable at lower temperatures, 
while at higher temperatures the stoichiometric (diamag- 
netic) structures become more stable. The temperature 
at which this transition occurs is a rapidly varying func- 
tion of po2- We have thus identified a class of sys- 
tems that undergo an unusual paramagnetic-diamagnetic 
transition induced by T and p of the reactive atmo- 
sphere, where the change in magnetic behavior reflects 
the change in composition. All the cluster sizes we 
have studied present this transition, but only for sizes 
3 < M < 8 the transition is predicted to happen in the 
pressure range 1 Pa < po^ < 1 atm. The transition be- 
tween paramagnetic and diamagnetic behavior is smooth 
when environmental conditions are changed smoothly, 
because different stoichiometrics always coexist within 
few fceT. In fact, the sand-colored areas in Fig. [3] are the 
regions where the free energies of the competing compo- 
sitions/structures are within 2k^T. 

In conclusion, we have presented a theoretical frame- 
work for predicting structure and stoichiometry of stable 
and metastable clusters in thermodynamic equilibrium 
with a gas atmosphere. A massively-parallel cascade ge- 
netic algorithm, i.e., an efficient scan of the potential 
energy of the clusters with various compositions, is com- 
bined with ah initio atomistic thermodynamics, a tool for 
evaluating the relative free energy of structures by know- 
ing their electronic energy, vibrational frequencies, and 
structural parameters for the evaluation of translational 
and rotational entropy. The methodology has been ap- 
plied to Mg clusters in an oxygen atmosphere. We have 
shown that small alkaline earth metal clusters form ther- 
modynamically stable non-stoichiometric "nano-oxides" , 
which have been overlooked so far. They are also pre- 
dicted to have peculiar (p,T)-dependent, magnetic prop- 
erties. 
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